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Abstract 

A  novel  algorithm  to  accurately  simulate  realistic  plasmas  with  constituent  species  at 
multiple  temperatures  using  an  extended  MHD  model  was  developed.  The  algorithm  was 
based  on  a  Roe-type  approximate  Riemann  solver.  The  algorithm  was  implemented  in  a  code 
to  model  the  time- dependent,  three-dimensional,  arbitrary-geometry  MHD  model  which  in¬ 
cludes  viscous  and  resistive  effects.  The  code  was  extended  to  include  thermal  diffusion  for 
the  constituent  temperatures  (neutrals,  ions,  and  electrons).  A  time-dependent  ionization 
model  was  added  which  self-consistently  calculates  the  ionization  fraction  of  the  fluid.  En¬ 
ergy  loss  mechanisms  were  added  for  the  constituent  fluid  components.  The  algorithm  was 
implemented  on  parallel  supercomputers  to  allow  the  detailed  modeling  of  realistic  plasmas 
in  complex  three-dimensional  geometries. 
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1  Executive  Summary 

The  primary  objective  of  this  project  is  to  develop  a  novel  algorithm  to  accurately  simulate 
realistic  plasmas  with  constituent  species  at  multiple  temperatures  using  an  extended  MHD 
model.  A  viable  time-dependent,  three-dimensional  MHD  code  will  provide  a  valuable  tool 
for  the  design  and  testing  of  plasma  related  technologies  that  are  important  to  the  Air  Force 
and  industry,  such  as  portable  pulsed  power,  high  power  microwave  devices,  hypersonic 
drag  reduction,  advanced  plasma  thrusters  for  space  propulsion,  nuclear  weapons  effects 
simulations,  radiation  production  for  counter  proliferation,  and  fusion  for  power  generation. 
Implementing  the  algorithm  on  parallel  supercomputers  will  allow  the  detailed  modeling  of 
realistic  plasmas  in  complex  three-dimensional  geometries. 

Current  MHD  codes  are  limited  to  simulations  of  short  time  scale  phenomena  because  of 
explicit  time  step  stability  limitations  and  equation  decoupling.  We  propose  developing  an 
implicit  algorithm  with  the  capability  to  simulate  physics  of  any  length  time  scale  because 
the  time  step  is  chosen  by  the  user  to  match  the  physics  of  interest.  This  algorithm  has  the 
additional  advantage  that  the  equations  are  solved  in  a  fully  coupled  manner.  The  plasma 
is  assumed  to  be  composed  of  a  neutral  fluid,  ion  fluid,  and  electron  fluid.  Each  fluid  has  an 
associated  temperature  and  can  exchange  energy  to  the  other  fluids  by  ionization  and  other 
collisional  processes.  The  plasma  is  allowed  to  have  a  variable  degree  of  ionization,  from  a 
fully  ionized  plasma  to  a  completely  neutral  gas.  The  algorithm  will  be  implemented  using 
arbitrary  finite  volumes  so  it  can  model  realistic  three-dimensional  geometries. 

To  speed  development  of  this  effort  an  existing  MHD  code  was  used.  WARPS  (Wash¬ 
ington  Approximate  Riemann  Plasma  code  for  3-d  domains)  is  a  time- dependent,  three- 
dimensional,  arbitrary- geometry  MHD  algorithm  with  viscous  and  resistive  effects.  The 
code  was  extended  to  include  thermal  diffusion  for  the  constituent  temperatures  (neutrals, 
ions,  and  electrons).  A  time-dependent  ionization  model  was  added  which  self-consistently 
calculates  the  ionization  fraction  of  the  fluid.  Energy  loss  mechanisms  were  added  for  the 
constituent  fluid  components.  These  features  were  benchmarked  against  analytical  results. 
Future  plans  include  the  addition  of  other  energy  loss  mechanisms  and  the  ability  to  transfer 
energy  between  constituent  fluid  components.  The  new  physics  algorithms  will  then  be  in¬ 
corporated  into  an  implicit  formulation  and  solved  using  a  domain  decomposition  technique 
for  parallel  computation. 

The  implicit  formulation  has  been  developed  for  the  resistive  and  viscous  MHD  model. 
The  culmination  of  this  research  effort  produced  the  Ph.  D.  dissertation  of  B.  Udrea.[l] 
The  algorithm  has  been  cast  using  finite  volumes  which  significantly  reduces  transverse  flux 
errors.  An  important  result  of  this  work  is  the  development  of  a  2- level  nested  iteration 
technique  which  accurately  solves  the  MHD  equations  with  typical  Courant  numbers  of  100. 
The  residual  of  the  error  is  driven  to  machine  accuracy  for  all  cases  investigated. 

As  a  result  of  this  project  several  professional  collaborations  now  exist  between  the 
Department  of  Aeronautics  and  Astronautics  at  the  University  of  Washington  and  the  Air 
Force  Research  Laboratory,  Lawrence  Livermore  National  Laboratory,  the  University  of 
Michigan,  the  University  of  Colorado,  Stanford  University,  and  other  departments  at  the 
University  of  Washington.  The  work  from  this  project  has  been  presented  at  international 
conferences  and  published  in  a  refereed  journal. 


2  Project  Description 

Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force,  some  of 
which  have  dual-use  potential.  These  applications  include  portable  pulsed  power  systems, 
high  power  microwave  devices,  drag  reduction  for  hypersonic  vehicles,  advanced  plasma 
thrusters  for  space  propulsion,  nuclear  weapons  effects  simulations,  radiation  production 
for  counter  proliferation,  and  fusion  for  power  generation.  Several  of  these  applications  are 
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specifically  mentioned  in  the  New  World  Vistas  Report  from  the  USAF  Scientific  Advisory 
Board.  [2]  In  general,  plasmas  fall  into  a  density  regime  where  they  exhibit  both  collective 
(fluid)  behavior  and  individual  (particle)  behavior.  Many  plasmas  of  interest  can  be  modeled 
by  treating  the  plasma  like  a  conducting  fluid  and  assigning  macroscopic  parameters  that 
accurately  describe  its  particle-like  interactions.  Magnetohydrodynamic  models  the  plasma 
in  this  manner. 

2.1  Research  Objectives 

The  objectives  of  the  project  are  to: 

•  Develop  an  implicit,  conservative  multi- temperature  algorithm  for  three-dimensional 
non-ideal  MHD  simulations  for  time-dependent  and  steady  state  variably  ionized  plas¬ 
mas; 

•  Validate  the  code  with  analytical  and  experimental  data;  and 

•  Apply  the  code  to  analyze  plasma  related  topics  at  the  Air  Force  Research  Laboratories 
[the  magnetic  flux  compression  generator  (MCG)  experiments  and  the  liner  implosion 
system  (WFX)  [3]  at  Kirtland  AFB ,  the  plasma  thruster  work  at  Edwards  AFB,  and  the 
hypersonic  drag  reduction  research  at  Wright-Patterson  AFB]  and  at  the  University  of 
Washington  [Z-Pinch  experiment  (ZaP)  [4]  and  Helicity  Injected  Tokamak  (HIT)  [5]]. 


2.2  Technical  Description 

The  three-dimensional,  extended  MHD  plasma  model  is  a  set  of  mixed  hyperbolic  and 
parabolic  equations.  The  Navier-Stokes  equations  are  also  of  this  type.  This  project  applies 
some  advances  that  have  been  made  in  implicit  algorithms  for  the  Navier-Stokes  equations 
to  the  MHD  equations.  These  implicit  algorithms  solve  the  equation  set  in  a  fully  cou¬ 
pled  manner,  which  generates  better  accuracy  than  the  current  methods  used  for  MHD 
simulations. 

When  expressed  in  conservative,  non-dimensional  form,  the  MHD  model  is  described  by 
the  following  equation  set. 
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pv 
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The  variables  are  density  (p),  velocity  (v),  magnetic  induction  (B),  ionization  fraction  (/*), 
pressure  (p),  and  energy  density  (e).  The  energy  density  is 
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v  v  B  B 

“b  P  n  +  — 


(2) 


where  7  =  cp/cv  is  the  ratio  of  the  specific  heats.  The  pressure  is  the  total  material  pressure, 
which  is  the  sum  of  the  partial  pressures  from  the  neutral,  ion,  and  electron  fluids. 


p  =  unkTn  -b  riikTi  +  nekTe 


(3) 


where  k  is  the  Boltzmann  constant,  nn  is  the  neutral  number  density,  and  Tn  is  the  neu¬ 
tral  temperature.  The  remaining  variables  are  the  ion  and  electron  number  densities  and 
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temperatures. 


P  =  (n„  +  m)Mi 


The  number  densities  are  determined  from  a  time-dependent  ionization  model 
d/n 

— —  =  7le  [<  (TV  >ion  —  ^  CTV  >recomb  ^i]  i  (^) 

at 

where  <  av  >ion  is  the  ionization  rate  parameter  and  <  av  >recomb  is  the  recombination 
rate  parameter. [6]  The  multiple  temperatures  evolve  independently  based  on  the  appropriate 
components  of  the  energy  equation  and  energy  transfer  between  species. 

The  right  hand  side  of  eqn(l)  contains  the  non-ideal  effects.  These  effects  include  vis¬ 
cosity,  resistivity,  Hall  currents,  diamagnetic  currents,  thermal  conduction,  and  radiation 
cooling.  The  non-ideal  terms  are  defined  by 

pVvisc  S  R^Al^  '  T  ^ 


— 


RmAl 


V  x  ^(VxB)) 


r>  __  WceTe  ^-7  w  X  B)  X  B 
Bhm  =  ~~RmAl  X  [  ~p 


Prad  s  -CTaiZ£ff{pfi)2TlJ2  (14) 

Mi  is  the  ion  mass,  ojcere  is  the  Hall  parameter,  Crad  is  the  Bremsstrahlung  radiation 
constant,  and  Zeff  is  the  effective  ionization  level  due  to  plasma  impurities. 

The  non-dimensional  tensors  are  the  stress  tensor  (f),  the  electrical  resistivity  (rj)>  and 
the  thermal  conductivity  (k),  and  I  is  the  identity  matrix.  The  non-dimensional  numbers 
are  defined  as  follows: 

Alfven  Number  : 

Reynolds  Number  : 

Magnetic  Reynolds  Number 
Peclet  Number  : 


Al  =  VA/V 
Re  =  LVjv 
Rm  =  fioLV/rj 
Pe  =  LV/k 
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The  characteristic  variables  are  length  (L),  velocity  ( V ),  Alfven  speed  ( Va  =  B/^/p0p), 
kinematic  viscosity  (i/),  electrical  resistivity  (77),  and  thermal  diffusivity  (k  =  k/pcp).  is 
the  permeability  of  free  space  (47r  x  10"7). 

For  convenience,  the  MHD  equation  set  [eqn(l)]  is  rewritten  in  the  following  compact 
form 

^  +  V-Th=QNon-ideal  (16) 

at 

where  Q  is  the  vector  of  conservative  variables,  is  the  tensor  of  hyperbolic  fluxes,  and 
Q Non-ideal  contains  the  non-ideal  (mostly  parabolic)  terms.  The  forms  of  these  vectors  and 
tensors  can  be  seen  from  the  previous  equations.  The  hyperbolic  fluxes  are  associated  with 
wave-like  motion,  and  the  parabolic  fluxes  are  associated  with  diffusion-like  motion. 


2.3  Multiple  Temperature  Evolution 

The  temperatures  of  the  multiple  species  (neutrals,  ions,  electrons)  evolve  independently 
based  on  the  appropriate  components  of  the  energy  equation  and  energy  transfer  between 
species.  The  temperatures  will  then  be  consistent  with  energy  conservation. 

The  temperature  rise  in  each  specie  will  depend  on  the  heating  mechanism  and  the 
density  fraction  of  the  specie.  We  define  the  density  fractions  as 
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and 
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(18) 


fe  =  fi- 


(19) 


For  multiply  charged  species,  the  last  definition  would  be  modified  to  fe  —  Zfc.  Viscous 
drag  heats  the  neutral  gas  and  the  ion  fluid,  but  does  not  affect  the  electrons.  Therefore, 
the  energy  rise  due  to  viscous  effects  is  attributed  to  the  neutral  and  ion  temperatures. 

.  =  (7-1  )fnPvisc  (20) 
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Resistivity  directly  heats  only  the  electrons. 

=  (7  -  1  )Pres 
Radiation  is  emitted  by  the  electrons  as  they  cool. 

=  (7  —  1 )  Prad 
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Each  species  has  its  own  thermal  conduction  component  as  is  evident  from  eqn(13). 

=  (7  -  1  )P*mJn  (24) 
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dpi 


dt 


cond 
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cond 

The  remaining  total  material  pressure  rise  is  due  to  adiabatic  compression  which  affects  all 
species  in  proportion  to  their  densities. 

Energy  may  also  transfer  between  species  due  to  collisions.  However,  total  energy  is 
conserved  and  the  total  material  pressure  will  not  be  affected.  The  interspecies  energy 
transfer  is  modeled  to  preserve  the  total  material  pressure. 

P  =  Pn  +  Pi  4-  Pe  =  Pn  +  Pi  +  Pe  (27) 


dpe 

dt 


2.4  Conservative  Algorithm 

Because  of  the  natural  differences  between  hyperbolic  and  parabolic  equations,  the  methods 
for  solving  them  are  very  different.  Since  the  MHD  equations  are  of  mixed  type  the  hyper¬ 
bolic  and  parabolic  terms  must  be  handled  differently.  The  hyperbolic  fluxes  are  differenced 
by  applying  an  implicit,  approximate  Riemann  algorithm  that  properly  accounts  for  their 
wave-like  behavior.  The  parabolic  terms  are  discretized  by  applying  explicit  central  differ¬ 
encing.  The  remaining  non-ideal  terms  which  correspond  to  the  Hall  effect  are  solved  using 
a  semi-implicit  method.  [7] 

The  design  of  the  algorithm  is  driven  by  the  conservative  numerical  techniques  that 
must  be  used  for  the  hyperbolic  terms.  Therefore,  we  begin  by  considering  the  ideal  MHD 
equations,  which  are  obtained  from  eqn(16)  by  setting  all  the  non-ideal  terms  (Q  Non-ideal) 
to  zero.  Note  that  ideal  MHD  refers  to  an  ideal  plasma  —  one  that  is  inviscid  and  non- 
resistive  and  neglects  thermal  conduction  and  finite  Larmor  radius  (FLR)  effects. 

The  ideal  MHD  equations  are 

+  =  +  -Q  =  0,  (28) 

dt  dt 


where  A  is  the  Jacobian  of  the  hyperbolic  flux  tensor. 


Here,  Q  is  the  vector  of  conserved  variables. 

rp 

Q  =  (p}pvXlpvy,pvz,BXiBy,Bz,e)  . 


(29) 


(30) 


This  is  a  set  of  hyperbolic  equations  and  thus  Ax  has  a  complete  set  of  real  eigenvalues 
given  by 

A  =  ( Vx,  0,  Vx  i  Vfast:  Vx  lb  Vslow>  V x  dz  Vax)  >  (31) 

where  Vfast  and  Vsiow  are  the  fast  and  slow  magnetosonic  speeds,  and  Vax  Is  the  Alfven 
speed  based  on  the  x  component  of  the  magnetic  field.  These  can  be  expressed  as 


V2  =  - 

*  fast  2 


C?  +  Vl  + 


'M+vF2 


(32) 


V,  =  - 

v  slow  2 


+  Vl  -  V/(C2  +  ^)2-4civX 


(33) 
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(34) 


Vj  - 


/  Ax 


flop' 


Here,  cs  is  the  ion  sound  speed,  which  for  a  perfect  gas  is 


(35) 


We  make  special  note  of  the  zero  eigenvalue,  A2  in  this  case.  The  zero  eigenvalue  only 
appears  in  multiple  dimensions  and  is  caused  by  the  perpendicular  nature  of  the  j  x  B  force. 
Powell  et  at,  [8]  recently  solved  this  zero  eigenvalue  problem  by  introducing  a  source  term 
that  is  proportional  the  divergence  of  the  magnetic  field.  The  eigenvalue  becomes  finite, 
A2  =  vx  in  this  case.  We  have  implemented  this  modification  and  it  is  discussed  in  a  later 
section. 

For  hyperbolic  equations  information  propagates  along  characteristics  which  travel  at 
wave  speeds  given  by  the  eigenvalues.  Most  modern  numerical  techniques  for  solving  hy¬ 
perbolic  equations  are  based  upon  the  idea  of  splitting  the  fluxes  into  components  due  to 
left-  and  right-running  waves.  Then  each  part  of  the  flux  can  be  differenced  in  an  upwind 
manner,  which  greatly  reduces  numerical  oscillations  and  stabilizes  the  solutions. 

It  is  well  known  that  if  a  hyperbolic  equation  is  solved  with  an  explicit  scheme,  then  the 
allowable  time  step  to  maintain  numerical  stability  is  given  by  the  CFL  (Courant-Friedrichs- 
Lewy)  condition,  which  in  the  case  of  the  ID  MHD  equations  is 


At  < 


Ax 

\vx  -f-  V/astj 


(36) 


For  the  high  magnetic  fields  and  low  densities  common  in  many  plasma  experiments,  the 
fast  magnetosonic  speed  is  quite  high,  and  thus  the  time  step  is  prohibitively  small.  We  are 
often  interested  in  only  modeling  the  physics  that  occurs  slower  than  Alfven  time  scales. 
For  example,  it  can  be  shown  that  resistive  tearing  modes,  which  are  important  in  studying 
fusion  plasmas,  evolve  on  a  time  scale  that  is  given  by [9] 

Ttearing  OC  T^T^5  =  (Luf/S  TA.  (37) 

ta  is  the  Alfven  time,  rv  is  the  resistive  diffusion  time,  and  Lu  is  the  Lundquist  number, 
which  is  given  by 


Lu=^L  =  RmAl. 
rA 


(38) 


If  Lu  is  106,  which  is  typical  for  laboratory  plasmas  in  fusion  applications,  the  resistive 
tearing  time  is  approximately  4000  times  larger  than  the  Alfven  time.  By  treating  the 
hyperbolic  fluxes  implicitly  in  time,  the  stability  restriction  on  the  time  step  is  removed, 
and  the  solution  can  be  advanced  at  the  larger  resistive  tearing  time  step.  This  is  our 
motivation  for  proposing  an  implicit  scheme. 


2.5  Implicit  Formulation 


For  clarity,  the  algorithm  for  the  two-dimensional  ideal  MHD  equations  is  presented.  The 
extension  to  three  dimensions  is  straight  forward.  The  two-dimensional  equation  is 


dQ  OF  8G  = 

dt  dx  +  dy 

Eqn(39)  was  discretized  using  first  order  Euler  time  differencing  to  get 


(39) 


(Q?P-Q?j) 


At 


-Rij  (Qn+1)  =  -R% 


(40) 
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where  R  is 


Rij  =  Fi+\j 2,j  “  -^i-l/2j  +  Gi,j  +  l/2  “  Gij_i/2.  (41) 

Note  that  in  this  equation  and  all  that  follow  the  grid  metric  terms  (cell  areas  and  volumes) 
were  omitted  for  clarity.  Linearizing  R  as  follows: 


(42) 

where  dR/dQ  has  been  defined  as  the  differenced  flux  Jacobians. 


dR 

dQ 


dF 

dQ 


*+1/2, j 


(43) 


where  dF/dQ  is  the  flux  Jacobians  of  the  x  flux.  The  flux  Jacobians  can  be  calculated 
analytically  with  the  assumption  that  the  solution  values  do  not  change  rapidly,  or  the  flux 
Jacobians  can  be  calculated  numerically.  Analytical  calculations  were  used  previously  [10] 
and  produced  adequate  results.  Numerical  calculation  was  investigated  for  the  current 
project.  Two  methods  for  numerical  calculation  were  investigated.  First,  a  limit  formulation 
similar  to  the  definition  of  a  differential  was  used. 


8F  =  F(Q  +  t)-F{Q)  + 
dQ  e 


(44) 


for  small  e  which  gave  first  order  accuracy  in  e.  The  flux  Jacobians  were  also  calculated 
using  complex  numbers.  The  flux  Jacobians  were  expanded  about  Q  in  a  Taylor  series. 


F(Q  +  ih)  =  F(Q)  +  ih^~ 


h2  d2F  .  h3  d3F 
T do2-1  6  dQ3 


+  ... 


(45) 


The  expression  was  rearranged  to  solve  for  the  flux  Jacobian. 


dF 

dQ 


F(Q  +  ih)' 
h 


+  0{h2) 


(46) 


which  gave  second  order  accuracy  in  h.  Additionally,  the  complex  formulation  required 
only  a  single  evaluation  of  the  flux  Jacobian  (though  using  complex  math)  compared  to  two 
evaluations  for  the  limit  formulation. 

Substituting  the  expression  for  back  into  eqn(40)  and  rearranging,  gave 


(47) 


Here  A Q  has  been  defined  as 

A  Qn  =  Ql+1-Q^. 


(48) 


The  left  hand  side  of  the  eqn(47)  is  an  implicit  operator  operating  on  A Q.  It  is  a  large 
banded  block  matrix.  In  three  dimensions,  it  is  an  (7max  ^  Jmax  X  i^max)  by  (I max  ^  Jmax  ^ 
Kmax)  matrix,  where  Imax  is  the  number  of  cells  in  the  x  direction,  etc.  It  is  quite  costly  to 
invert  a  matrix  of  this  size  directly.  For  this  project  a  more  efficient  iterative  method  was 
chosen.  When  solved  iteratively  ,  eqn(47)  can  lose  time  accuracy.  To  recover  time  accuracy 
the  time  derivative  of  Q  was  added  as  a  source  term  to  the  right  hand  side  of  the  equation. 
The  modified  equation  became 


© 


nn+l  cn  +  l 

^ ij  ^ ij 


(49) 
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where 

sr/,  =  5^(3«”/1-4orJ+«r,)“§-  (“) 

The  r  in  eqn(49)  is  an  iteration  variable, called  pseudo  time.  At  each  physical  time  step, 
eqn(49)  is  solved  iteratively  until  the  left  hand  side  vanishes.  When  the  solution  converges, 
our  original  equation 


dQ 

dt 


=  -R 


(51) 


is  solved.  This  technique  is  known  as  dual  time-stepping.  [11]  Note  that  in  eqn(50)  a  more 
accurate  time  derivative  can  be  used  at  the  expense  of  the  additional  memory  needed  to 
store  the  Q  vectors  from  previous  time  steps. 

One  advantage  of  the  strategy  outlined  above  is  that  the  implicit  operator  and  the  right 
hand  side  in  eqn(47)  are  decoupled.  The  structure  of  the  matrix  no  longer  depends  on  the 
details  of  the  discretization  of  the  right  hand  side  fluxes.  The  iteration  equation  [eqn(47)] 


31  1 1 

SdR\m' 

3Ar  -f-  2A t 

\dQJij 

A QTj  =  ~  (3 Q%  -  4 Q?j  +  Qif1)  -  R?i 


(52) 


has  the  standard  form 


Ax  =  b. 


(53) 


Previously  the  LU-SGS  method  (lower-upper  symmetric  Gauss-Seidel)[12]  was  used  to 
invert  the  implicit  operator  A.  [13,  10]  The  LU-SGS  method  required  a  modification  of  the 
implicit  operator  through  an  approximate  factorization  procedure  which  reduced  the  accu¬ 
racy  of  the  operator  and  led  to  poor  convergence. 

The  current  project  used  a  symmetric  Gauss-Seidel  method  which  does  not  rely  on 
approximate  factorization.  The  SGS  method  was  used  to  iteratively  invert  the  implicit 
operator  and  the  approximate  Riemann  solver  that  is  used  to  form  the  right  hand  side 
fluxes. 

The  implicit  operator  matrix  was  decomposed  into  lower,  upper,  and  diagonal  matrices 
and  written  as 


A  -  L  -f  U  4-  D  (54) 

Each  iteration  of  the  symmetric  Gauss-Seidel  method  performs  two  sweeps  of  the  domain 
—  a  forward  sweep  followed  by  a  backward  sweep. 


(L  +  D)x(2i~'1)  +  Ux(2i~2)  =  b 

(55) 

(U  +  D)x(2°  +  Lx(21_1)  =  b 

(56) 

where  Z  —  1,2,3,...  is  the  iteration  index. 

For  reasonably  large  values  of  Re  and  Rm  (easily  in  the  range  of  interest  for  most 
applications),  the  parabolic  terms  can  be  differenced  explicitly  without  constraining  the 
allowable  time  step.  The  right  hand  side  of  eqn(47)  was  modified  by  adding  the  central 
differenced  parabolic  terms. 
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2.6  Approximate  Riemann  Solver 

The  fluxes  on  the  right  hand  side  of  eqn(47)  were  discretized  using  a  Roe-type  approximate 
Riemann  solver.  [14]  In  this  method  the  overall  solution  was  built  upon  the  solutions  to  the 
Riemann  problem  defined  by  the  discontinuous  jump  in  the  solution  between  each  pair  of 
cells.  The  numerical  flux  for  a  first-order  accurate  (in  space)  Roe- type  solver  was  written 
in  symmetric  form  as 

Fi+ 1/2  =  \  (Fi+ !  +Fi)-\Y/lk  (Qi+ 1  -  Qi)  lAfc|  rk  (57) 

where  r*,  is  the  kth  right  eigenvector,  A*,  is  the  absolute  value  of  the  kth  eigenvalue,  and  Ik 
is  the  kth  left  eigenvector.  The  values  at  the  cell  interface  (2  4-1/2)  were  obtained  by  either 
a  simple  average  or,  more  accurately,  a  Roe-average  of  the  neighboring  cells.  Determining  a 
Roe-average  on  an  arbitrary  computational  grid  involved  transforming  the  vector  quantities 
to  a  coordinate  system  that  is  orthogonal  to  the  local  cell  interface.  Then  the  flux  calculated 
as  above  will  be  normal  to  the  cell  interface  which  is  the  desired  orientation  for  applying 
the  divergence  theorem  in  a  finite  volume  method. 

These  first  order  accurate  upwind  fluxes  are  used  in  the  vicinity  of  sharp  discontinuities 
in  order  to  suppress  oscillations  in  the  solution.  Globally  second  order  accurate  solutions 
were  achieved  by  using  a  flux  limiter  that  modifies  the  first  order  flux  so  that  it  uses  second 
order  central  differencing  in  smooth  portions  of  the  flow.  The  minmod  limiter  was  used. [15] 

Once  the  eigenvalues  and  eigenvectors  were  obtained  and  properly  normalized  to  avoid 
singularities,  it  was  relatively  straight-forward  to  apply  this  scheme  to  the  one-dimensional 
ideal  MHD  equations.  [16,  17]  Unlike  for  the  Euler  equations,  the  extension  to  more  than 
one  dimension  was  non-trivial.  In  more  than  one  dimension,  the  Q  vector  must  include  Bx 
in  addition  to  the  other  magnetic  field  components.  (For  the  one- dimensional  case  Bx  is 
constant  by  virtue  of  V  *  B  =  0).  Since  the  j  x  B  force  acts  perpendicularly  to  the  directions 
of  j  and  B,  the  F  flux  vector  has  a  zero  term  corresponding  to  Bx.  Thus,  the  Jacobian 
matrix  of  F  is  singular  and  has  a  zero  eigenvalue.  A  complete  set  of  physically  meaningful 
eigenvectors  no  longer  exists.  Physically,  one  would  expect  information  to  travel  either  at 
the  fluid  velocity  or  at  the  fluid  velocity  plus  or  minus  the  wave  speeds.  Simply  dropping 
Bx  from  the  equation  set  is  not  a  viable  option,  because  Bx  needs  to  change  in  order  to 
maintain  V-B  =  0.  Powell  et  ai,[S]  recently  solved  this  problem  by  modifying  the  Jacobian 
in  such  a  way  as  to  change  the  zero  eigenvalue  to  vx  (keeping  the  others  unchanged),  and 
then  adding  in  a  source  term  that  exactly  canceled  out  the  terms  introduced  by  the  modified 
Jacobian. 

The  source  term  is 


Sdi- 


P 

B 

v 

v  B 


V-B 


(58) 


It  is  proportional  to  the  divergence  of  B  and  thus  very  small. 


2.7  Finite  Volume  Grid  and  Parallel  Implementation 

Since  the  algorithm  being  developed  will  be  used  for  real  devices,  it  must  be  capable  of 
modeling  arbitrary,  three-dimensional  geometries.  Therefore,  multi-block,  finite  volume 
grids  were  used.  A  typical  cell  is  shown  in  Figure  1.  The  computational  domain  is  divided 
into  blocks  which  are  then  spanned  by  body-fitting  finite  volume  cells.  See  Figure  2  for  a 
possible  grid. 
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Figure  1:  Schematic  of  the  arbitrary  shaped  three-dimensional  finite  volume  cell  used  by  the 
algorithm. 


Figure  2:  Schematic  of  the  arbitrary  shaped  three-dimensional  finite  volume  cell  used  by  the 
algorithm. 


As  discussed  in  the  previous  section,  the  formulation  of  the  approximate  Riemann  solver 
which  we  have  developed  generates  hyperbolic  fluxes  oriented  normal  to  the  grid  cell  in¬ 
terfaces.  Application  of  the  divergence  theorem  is  then  a  simple  operation.  The  parabolic 
fluxes  are  also  calculated  to  be  normal  to  the  grid  cell  interfaces.  To  accomplish  this  ori¬ 
entation,  a  set  of  nested  control  volumes  were  used  and  the  appropriate  vector  operations 
within  these  volumes  were  applied. 

The  block  structure  of  the  grid  provided  a  natural  domain  decomposition  for  the  parallel 
implementation.  The  integral  form  of  a  general  conservation  law  was  expressed  as 


|  J  dV  Q  +  dS  ■  F(Q) 
n  s 


J  dV  S(Q), 

n 


(59) 


where  ft  is  the  domain  and  £  is  the  boundary  of  ft.  Q  is  the  vector  of  conserved  variables, 
F{Q)  is  the  flux  of  the  conserved  variables,  and  S{Q)  is  the  vector  of  source  terms.  By 
splitting  the  domain  ft  into  p  subdomains  such  that 


v 

n  =  \Jnt,  (60) 

t=l 


eqn(59)  was  replaced  with  a  set  of  p  conservation  equations  applied  on  the  subdomains  ft*. 


jt  J  dVQ  +  j>dS-F{Q) 

fti  Si 


=  JdV  S(Q ). 

Ui 


i  =  1,2 


(61) 


Each  of  these  discretized  equations  is  solved  by  a  single  processor  using  the  boundary  values 
copied  from  neighboring  subdomains. 

To  ensure  a  portable  code  a  message  passing  system  commonly  available  on  parallel 
supercomputers  and  on  workstation  clusters  was  used.  This  system  was  the  Message  Pass¬ 
ing  Interface  (MPI)  [18],  which  was  adopted  as  a  standard  in  May  1994  by  industry  and 
academia.  Hardware  and  software  vendors’  implementation  of  MPI  provides  parallel  pro¬ 
gram  developers  with  a  consistent  set  of  subroutines  callable  from  FORTRAN90  and  C.  In 
this  project  the  basic  point-to-point  communications  subroutines  and  global  communica¬ 
tions  subroutines  were  used.  The  point-to-point  communication  subroutines  were  used  for 
the  domain  decomposition  and  boundary  exchange  while  the  global  communication  sub¬ 
routines  were  used  for  convergence  checking.  All  message  passing  systems  (PVM,  MPL) 
support  point-to-point  and  global  communications  subroutines  so  that  by  using  only  the 
basic  set  portability  to  systems  not  supporting  MPI  was  simplified. 


3  Project  Implementation  and  Results 

3.1  Finite  Volume  Improvements 

To  improve  the  codes  ability  to  handle  highly  distorted  grids,  finite  volume  grids  were 
implemented.  The  finite  volume  implementation  greatly  reduced  the  anomalous  momentum 
leakage  into  orthogonal  directions  when  the  grid  was  distorted. 

Figure  3  shows  a  shock  tube  test  problem.  The  simulation  was  performed  in  three- 
dimensions,  but  should  remain  one-dimensional.  The  figure  shows  the  three  gas  dynamic 
waves  in  the  density  plot.  The  transverse  velocity  should  be  zero.  A  finite  amount  of 
momentum  leakage  was  generated  by  the  grid  metrics  in  the  finite  difference  generalized 
coordinate  formulation. 

The  reduction  of  the  anomalous  momentum  leakage  into  orthogonal  directions  can  be 
seen  in  Figure  4.  The  simulation  was  identical  to  that  shown  in  Figure  3  except  a  finite 
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Figure  3:  A  three-dimensional  simulation  of  a  one-dimensional  shock  tube  showing  the  density 
and  transverse  velocity.  The  three  gas  dynamic  waves  can  be  seen  in  the  density  plot.  The 
transverse  velocity  should  be  zero. 
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Figure  4:  A  three-dimensional  simulation  of  a  one- dimensional  shock  tube  showing  the  density 
and  transverse  velocity.  The  three  gas  dynamic  waves  can  be  seen  in  the  density  plot.  The 
transverse  velocity  should  be  zero. 

volume  implementation  was  used  instead  of  the  generalized  coordinate  formulation.  Fig¬ 
ure  4  shows  the  same  three  gas  dynamic  waves  in  the  density  plot.  With  the  finite  volume 
implementation,  the  transverse  velocity  is  zero  to  within  machine  accuracy. 

3.2  Implicit  Formulation  and  Numerical  Flux  Jacobian  Calcula¬ 
tions 

The  LU-SGS  method  (lower-upper  symmetric  Gauss-Seidel)[12]  was  used  previously  to  in¬ 
vert  the  implicit  operator  of  eqn(47).[13,  10]  The  LU-SGS  method  required  a  modification 
of  the  implicit  operator  through  an  approximate  factorization  procedure  which  reduced  the 
accuracy  of  the  operator  and  led  to  poor  convergence.  The  convergence  history  is  shown 
in  Figure  5  The  explanation  for  the  poor  convergence  was  inaccurate  approximation  of 
the  implicit  operator.  The  inaccuracy  developed  from  the  combination  of  the  approximate 
factorization  and  the  approximate  analytical  calculation  of  the  flux  Jacobians. 

The  standard  formulation  of  eqn(52)  allowed  the  use  of  standard  iterative  matrix  in¬ 
version  methods.  In  this  project  the  symmetric  Gauss-Seidel  method  was  used.  The  con¬ 
vergence  history  is  shown  in  Figure  6.  n  is  the  number  of  physical  time  iterations,  m  is 
the  number  of  pseudo  time  subiterations,  and  sgs  is  the  number  of  iterations  of  the  SGS 
method.  Unlike  the  LU-SGS  method,  there  is  no  coupling  between  the  SGS  iterations  and 
the  pseudo  time  iterations.  The  values  of  the  implicit  operator  A  and  the  inhomogeneity 
vector  b  are  updated  between  pseudo  time  iterations,  but  not  between  SGS  iterations. 
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Figure  6:  Convergence  history  using  the  SGS  method  to  invert  the  implicit  operator,  n  is  the 
number  of  physical  time  iterations,  m  is  the  number  of  pseudo  time  subiterations,  and  sgs  is 
the  number  of  iterations  of  the  SGS  method. 
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Figure  7:  Mach  3  flow  impinging  on  a  hemispherical  body.  The  upper  plot  shows  the  density  and 
streamlines.  The  lower  plot  shows  the  ionization  state  and  streamlines.  Notice  the  increased 
ionization  at  the  stagnation  point. 


The  formulation  of  eqn(52)  into  a  standard  form  required  numerical  calculation  of  the 
flux  Jacobians.  The  limit  calculation  given  by  eqn(44)  was  simple  to  implement  and  gave 
accurate  values  of  the  flux  Jacobian.  However,  it  was  sensitive  on  the  value  of  e  when  e 
was  increased  beyond  1  x  10-12.  Using  the  complex  formulation  given  by  eqn(46)  provided 
accurate  flux  Jacobian  calculations  with  much  less  sensitivity  on  h.  Typical  values  of  h  were 
1  x  1CT5. 

3.3  Time  Dependent  Ionization  and  Multiple  Temperature  Effects 

A  time-dependent  ionization  model  was  added  to  self- consistently  calculate  the  ionization 
fraction  of  the  fluid.  The  model  is  described  by  eqn(6).  The  ionization  rate  parameter 
<  c tv  >ion  and  the  recombination  rate  <  av  > recomb  were  calculated  using  empirical  formu¬ 
lations  given  in  Ref.  [6]. 

The  time- depen  dent  ionization  model  allowed  determination  of  the  ionization  fraction. 
For  uniform  flow  properties  the  time  dependence  is  exponential.  The  model  was  bench- 
marked  against  analytical  formulations  for  its  time-dependence.  The  steady-state  solution 
was  benchmarked  against  the  Saha  equation. 

A  more  interesting  test  was  constructed  to  have  a  Mach  3  flow  impinge  on  a  hemispherical 
body.  The  flow  was  initially  unionized,  and  ionized  upon  transition  through  the  shock  wave. 
The  results  are  shown  in  Figure  7. 

A  subset  of  the  multi-temperature  effects  have  been  added  to  the  code.  The  code  was 
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Figure  8:  (a)  Strip  decomposition  and  (b)  patch  decomposition  of  a  2-D  domain. 

extended  to  include  thermal  diffusion  for  the  constituent  temperatures  (neutrals,  ions,  and 
electrons).  Energy  loss  mechanisms  were  added  for  the  constituent  fluid  components.  The 
constituent  temperatures  can  evolve  independently  by  diffusion  as  described  in  eqn(13). 
The  heat  conduction  equation  was  used  as  a  benchmark  to  validate  the  incorporation  of  the 
multi-temperature  thermal  diffusion  model. 

A  dominate  energy  loss  mechanism  for  high  electron  temperature  plasmas  is  radiation. 
The  radiation  loss  term  due  to  Bremsstrahlung  has  been  included.  The  radiation  was 
validated  against  analytical  results. 

3.4  Parallel  Computer  Performance 

The  algorithm  was  parallelized  using  the  domain  decomposition  technique  (DDT).  The 
integral  form  of  a  general  conservation  law  was  given  by  eqn(59).  Domain  decomposition 
implementation  requires  boundary  (or  ghost)  cells  to  overlap  with  neighboring  domains  or 
blocks.  The  domain  decomposition  is  illustrated  for  two  dimensions  in  Figure  8. 

The  ghost  cells  are  used  as  boundary  conditions  to  the  real  cells  in  the  block.  A  con¬ 
sequence  of  domain  decomposition  is  the  more  blocks  that  are  used  the  more  ghost  cells 
that  are  necessary.  The  ghost  cell  data  lag  the  current  computation  by  a  single  iteration. 
Therefore,  an  increase  on  ghost  cells  generated  a  slower  convergence  rate.  Figure  9  shows 
the  slightly  slower  convergence  rate.  The  convergence  history  for  4  and  8  processors  overlay. 

The  parallel  performance  for  the  code  is  shown  in  Figures  10  and  11  for  the  code  operat¬ 
ing  in  explicit  and  implicit  mode.  The  grid  was  the  three-dimensional  grid  shown  in  Figure  2 
and  was  parallelized  using  domain  decomposition.  The  grid  was  scaled  with  the  number  of 
processors,  so  the  grid  size  per  processor  was  constant.  As  the  number  of  processors  was 
doubled,  the  number  of  grid  cells  was  also  doubled.  The  ideal  speedup  was  unity.  Note  that 
the  speedup  presented  is  “engineering”  speedup.  The  value  includes  not  only  the  inefficien¬ 
cies  associated  with  communication  between  the  processors  but  also  those  associated  with 
more  iterations  required  to  converge  the  solution.  The  “engineering”  speedup  is,  therefore, 
the  total  parallel  efficiency  to  obtain  the  same  quality  of  solution  on  a  parallel  computer. 

The  high  parallel  efficiency  was  obtained  by  overlapping  communication  with  compu¬ 
tation.  The  boundary  information  was  exchanged  while  the  core  cells  were  computed  and 
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Figure  9:  Convergence  history  using  the  SGS  method  to  invert  the  implicit  operator  on  a 
parallel  computer.  The  results  from  a  serial  computer  are  plotted  for  comparison,  n  is  the 
number  of  physical  time  iterations,  m  is  the  number  of  pseudo  time  sub  iter  at  ions,  and  sgs  is 
the  number  of  iterations  of  the  SGS  method. 
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Figure  10:  Parallel  speedup  for  a  three-dimensional  grid  using  domain  decomposition  on  a 
cluster  of  DEC  Alpha  workstations.  The  grid  is  scaled  with  the  number  of  processors,  so  the 
grid  size  per  processor  is  constant.  The  ideal  speedup  is  unity. 

before  the  boundary  cells  are  computed.  The  super-linear  speedup  for  the  explicit  operating 
mode  was  generated  by  slow  operation  on  a  single  processor.  This  effect  has  been  reported 
Michl  et  al  [19]  Figure  10  contains  the  performance  results  from  a  parallel  workstation  clus¬ 
ter  of  16  DEC  Alpha  workstations.  Figure  11  contains  the  performance  results  from  the 
IBM  SP2  at  the  Maui  High  Performance  Computing  Center.  The  results  on  both  computing 
platforms  were  similar. 

4  Professional  Interactions 

4.1  Project  Personnel 

The  personnel  who  have  been  directly  involved  in  this  project  are  listed  below. 

_ Name _ Position _ 

Uri  Shumlak  Assistant  Professor 

D.  Scott  Eberhardt  Associate  Professor 

Thomas  R.  Jarboe  Professor 

R.  Scott  Raber  Graduate  Student 

Bogdan  Udrea  Graduate  Student 

Ward  Vuillemot  Graduate  Student 

Graham  Schelle  Undergraduate  Student 

4.2  Collaborations 

4.2.1  Air  Force  Research  Laboratory 

Dr.  Robert  Peterkin,  Jr.  of  the  Electromagnetic  Sources  Division  of  the  Air  Force  Research 
Laboratory  at  Kirtland  AFB  on  three-dimensional  multigrid  algorithms  for  MACH3,  devel¬ 
opment  of  a  parallel  PIC  (particle  in  cell)  code  for  microwave  simulations,  and  stabilization 
of  the  Rayleigh-Taylor  instability  in  solid  liner  implosions  by  introducing  a  sheared  axial 
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No.  of  Processors 


Figure  11:  Parallel  speedup  for  a  three-dimensional  grid  using  domain  decomposition  on  IBM 
SP2  parallel  supercomputer.  The  grid  is  scaled  with  the  number  of  processors,  so  the  grid  size 
per  processor  is  constant.  The  ideal  speedup  is  unity. 


flow.  A  full  approximation  multigrid  method  was  installed  in  MACH3  to  calculate  thermal 
diffusion  and  magnetic  resistive  diffusion.  Knowledge  developed  in  the  area  of  relaxation 
schemes  was  implemented  into  the  ICEPIC  code  to  make  a  3-D  Poisson  solver.  The  solver 
was'  needed  to  determine  electric  field  concentration  on  a  high  power  microwave  source. 
The  collaboration  occured  in  person  during  July  and  August.  Several  phone  and. email 
discussions  took  place  throughout  the  year. 

4.2.2  Sandia  National  Laboratories 

Dr.  Norm  Roderick  of  the  Pulsed  Power  Sciences  Center  at  Sandia  National  Laboratories 
on  the  uses  of  sheared  axial  flows  to  stabilize  z-pinch  implosions.  This  is  an  ongoing  collab¬ 
oration  that  resulted  in  the  publication  listed  in  the  following  section. 

4.2.3  University  of  Washington 

Prof.  Scott  Eberhardt  of  the  Aeronautics  and  Astronautics  Department  and  Prof.  Randy 
LeVeque  of  the  Applied  Math  Department  on  approximate  Riemann  solvers  and  their  ap¬ 
plications  to  multidimensional  problems.  We  have  regular  discussions  on  a  weekly  basis. 

Prof.  Tom  Jarboe  of  the  Aeronautics  and  Astronautics  Department  on  the  higher  mode 
stability  of  spheromaks  and  on  the  effect  of  realistic  three-dimensional  geometries  on  sphero- 
mak  stability.  This  is  an  ongoing  collaboration  that  resulted  in  the  publications  listed  in 
the  following  section. 

4.3  Publications 

A  journal  article  describing  our  algorithm  has  been  published  in  the  Journal  of  Computa¬ 
tional  Physics.  The  title  is  “An  Implicit  Scheme  for  Nonideal  Magnetohydrodynamics”  by 
0.  S.  Jones,  U.  Shumlak,  and  D.  S.  Eberhardt. [13]  The  citation  is  Journal  of  Computational 
Physics  130,  231  (1997).  Another  journal  article  which  describes  the  use  of  sheared  flows  to 
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stabilize  the  Rayleigh-Taylor  instability  has  been  published  in  Physics  of  Plasmas.  This  is 
work  that  was  performed  with  collaboration  at  the  Air  Force  Research  Laboratory.  The  title 
is  “Mitigation  of  the  Rayleigh-Taylor  Instability  by  Sheared  Axial  Flows”  by  U.  Shumlak 
and  N.  F.  Roderick. [20]  The  citation  is  Physics  of  Plasmas  5,  2384  (1998). 

Two  papers  describing  the  higher  mode  stability  in  spheromak  plasmas  and  the  effect  of 
realistic  three-dimensional  geometries  on  spheromak  stability  have  also  been  published. [21, 
22]  The  citation  is  Physics  of  Plasmas  6,  4382  (1999).  The  second  paper  has  been  submitted 
for  publication  in  Physics  of  Plasmas. 

4.4  Presentations 

A  paper  was  presented  at  the  Forty  First  Annual  American  Physical  Society  Meeting  of  the 
Division  of  Plasma  Physics,  Seattle,  Washington,  November  1999.  The  title  was  “Plasma 
Effects  on  Hypersonic  Flows,”  by  W.  Vuillemot  and  U.  Shumlak; 

5  Conclusions 

The  successful  development  of  the  three-dimensional  advanced  implicit  algorithm  and  the 
implementation  of  time-dependent  ionization  and  multiple  temperature  effects  show  that 
this  project  is  reaching  its  objectives.  The  research  related  to  this  project  has  been  published 
in  refereed  journals  and  presented  at  international  conferences.  Valuable  collaborations  have 
been  formed  with  the  Air  Force  Research  Laboratory,  Sandia  National  Laboratory,  and  other 
universities. 

The  continuing  development  of  this  project  will  include  investigating  more  powerful  im¬ 
plicit  matrix  inversion  methods,  the  addition  of  other  energy  loss  mechanisms  and  the  ability 
to  transfer  energy  between  constituent  fluid  components.  The  new  physics  algorithms  will 
then  be  incorporated  into  an  implicit  formulation  and  solved  using  a  domain  decomposition 
technique  for  parallel  computation.  The  code  will  also  be  applied  to  plasma  experiments 
to  calibrate  the  code  and  gain  physical  insight  into  devices  that  are  important  to  the  Air 
Force. 
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